arXiv:1506.08646v2 [hep-ph] 6Jul2015 


Evolution of heavy quark distribution function in quark-gluon 
plasmas: using the Iterative Laplace Transform Method 


Sharareh Mehrabi Pari * * 1,2 , Kurosh Javidan T and Fatemeh Taghavi Shahri ^ 

department of Physics, Ferdowsi University of Mashhad, 91775-1436 Mashhad, Iran 
department of Physics and Technology, University of Bergen, 5007 Bergen, Norway 


Abstract 

The ’’Iterative Laplace Transform Method” is used to solve the Fokker-Planck 
equation for finding the time evolution of the heavy quarks distribution functions 
such as charm and bottom in quark gluon plasma. These solutions will lead us to 
calculation of nuclear suppression factor Raa- The results have good agreement with 
available experiment data from the PHENIX collaboration. 


I. Introduction 

Colliding heavy ions in Relativistic Heavylion Collider (RHIC) and Large hadron Col- 
lider(LHC) experiments generate temperatures many thousand times the temperature of 
the sun. These collisions are recreating in the laboratory conditions similar to those just 
after the big bang. Under these extreme conditions, protons and neutrons melt, freeing the 
quarks from their bonds with the gluons. This new phase of matter is called quark-gluon 
plasma. This highly excited state of matter displays properties similar to a nearly per¬ 
fect fluid, which has been successfully described by hydrodynamic models [DEI SIS! -This 
medium is probed by hard processes with scales ranging from a few GeV to several of TeV. 
The existence of such a phase and its properties are key issues in the theory of quantum 
chromodynamics (QCD). Basically, heavy partons, which are produced in early stage of 
heavy ion collisions, have relatively high density. Such heavy partons could induce a large 
amount of energy loss for hard patrons produced in the initial stage of the collision. So, 
the observed jet quenching in such collisions can indicate the formation of this strongly 
interacting dense matter [5], El El El EE [lOj [IT] . Charm and bottom quarks play crucial 
role in this endeavor because they are the witness to the entire space-time evolution of 
the system and due to their large masses, a memory of their interaction history may be 
preserved [12l 13, l4l fl5] . 

Their mass is significantly larger than the temperature of the plasma they are formed 
in, M 2> T c . They are produced in the early stage of heavy ion collisions and they do 
not dictate the bulk properties of the matter. In addition the heavy quark thermalization 
time in perturbative QCD is estimated to be of order of 10 — 15 fm/c and 25 — 30 fm/c 
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for charm and bottom quarks respectively [16] [T71 18l T9] . This time is larger than the 
thermalization time of the light quarks and gluons. 

The plasma reaches local equilibrium at an initial temperature T) and initial thermal¬ 
ization time Tj. It is assumed that its constituent particles are light quarks, anti-quarks and 
gluons while nonthermalized heavy quarks, which are formed in early stages of collision, 
are particles with Brownian motion that move through the expanding QGP background. 
The Fokker-Planck equation (FPE) provides an appropriate framework for studying the 
evolution of the heavy quark distribution functions [20, ,21]. 

Heavy quarks dissipate their energy during their propagation through QGP via two 
processes [22]: 

The first one is the collisional process, such as gQ —> gQ, qQ —> qQ, qQ —> qQ 

The second one is the radiative process, such as Q + q —> Q + q + g, Q + g — > 
Q + g + g [221241125]. 

These processes can be described well by the Fokker-Plank equation. 

In this paper our main task is to study the evolution of heavy quarks in quark gluon 
plasma by solving the Fokker-Plank equation with new method named ’’Iterative Laplace 
Transform Method” Hi [27]. This method gives numerical solutions in the form of con¬ 
vergent series with easily computable components. To verify our solutions, we calculate 
the nuclear suppression factor, Raa{pt), of heavy quarks using simulated data of heavy 
ion collisions [28], 29 , 30] . 

To do this we consider both the collisional and radiation parts of energy loss in FPE 
and finally we compare our results with experimental data from PHENIX [ST]. 

This paper is organized as follows: In the next section we review the Fokker-Planck 
equation and the collisional energy loss of heavy quarks as well as their energy loss due to 
radiation. In section III we introduce iterative Laplace transform method. In Section IV 
we study the time evolution of heavy quarks in QGP using the presented method. Section 
V is focused on calculating the nuclear suppression factor Raa(pt)- The last section is 
devoted to concluding remarks. 


II.The heavy quark momentum evolution in FPE dynamics 


The Quark Gluon Plasma (QGP) is formed at the time Tj after the impact. Time evolution 
of the momentum distribution of non-equilibrated heavy quarks due to the interaction 
with the equilibrated QGP during the time interval r* < r < thq can be calculated by 
the Fokker-Planck equation. Evolution of the expanding QGP is described by relativistic 
hydrodynamics [32l 33]. 

To derive the Fokker-planck equation, let us start from the Boltzmann transport equa¬ 
tion, which in covariant form is [34] : 


P lx d fl f(x,p) = c{f } 


( 1 ) 


where p^ 1 is the four-momentum of the heavy quarks (HQs), C{/} is the collision term and 
f(x,p ) is the momentum distribution function of the heavy quarks. For uniform plasma, 
/ will lose its dependence on x , so the Boltzmann equation becomes 


df _c{f} _ r df^ 
dt E [ dt >coU 


( 2 ) 


To simplify the non-linear integro-differential Boltzmann equation, we employ the Lan¬ 
dau approximation. This means that we consider soft scatterings with small momentum 
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transfer compared to the particle momentum p, and neglect the interaction of heavy quarks 
with each other. According to these assumptions the evolution of momentum distribution 
of heavy quarks is given by 


d_l 

dt 


d_ 

dpi 




(3) 


where’d, 21; (p) and Q5;j(p) are drag and diffusion coefficients respectively, defined as: 


21; = / d 3 kuj(p,k)ki 


(4) 


*ij = 


1 [ „ 

- / d 6 kLo(p,k)kikj. 


(5) 


Here, we considered the elastic scattering of the HQs with gluons, light quarks, and cor¬ 
responding anti-quarks. All these interactions contribute to determine the collision rate 
w{p,k) [35]. For an isotropic QGP one can write 21; = PiA(p) and 25;^ = D(p)5ij. In 
the simplest form, we assume that the drag and diffusion coefficients are slowly varying 
functions of momentum, so Equation ([3]) reduces to the Fokker-Planck equation: 


df{p,t ) 

at 


A (p)-^(pf(p,t)) + D(p) 


d_ 

dp 


f{p,t ) 


( 6 ) 


The drag coefficient is an important quantity carrying information about the dynamics 
of elastic collisions. It is expected that the drag coefficient should be determined by the 
properties of the bath and not by the nature of the HQ [36j. 

After heavy ion collision, the charm quarks are produced at a time scale of about 
0.07 fm/c. This time is about 0.02 fm/c for bottom quarks. These heavy quarks will 
propagate through the deconfined matter (QGP) while losing energy. 

In this work we use the drag coefficient as follows: 


A(p,t) 


IdE 
p dL 


(7) 


If coupling between HQs and background QGP is weak, then the diffusion coefficient can be 
calculated using the Einstein relation, D(p) = MTA(p), where M is the mass of the heavy 
quark and T is the thermal bath temperature [57 , 38]. The first perturbative estimation 
of the energy loss in a QGP was proposed by Bjorken for heavy quarks (39|, which is an 
important quantity. It may be noted that the jet quenching, drag force, damping rate 
of particles in the plasma and etc., are determined directly in terms of the energy loss. 
Furthermore, since the calculated energy loss is proportional to the heavy fermion mass 
in quark-gluon plasma, it is used as an important tool to identify particles and give more 
information about the properties of QGP. The collisional energy loss of a heavy quark with 
energy E and mass A 1 in a thermalized medium with temperature T has been calculated 
in several papers with different approaches [3D]. In our calculations, the collisional energy 
loss of HQs in QGP by considering ’’Hard and Soft Thermal Loops” are given as follows 

m- 
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where v is the HQ speed, ot s is the strong coupling constant, n, is the number of quark 
flavors in the medium, E and M are energy and mass of the propagated HQ respectively, 
m g = yj{ 1 + rif/6)g 2 T' 2 /3 is the thermal gluon mass, g = VTkcT s is the gauge coupling 
parameter and B{y) is a smooth velocity function, which can be taken approximately as 
0.7 [H]. 


The collisional energy loss of a HQ in a thermalized QGP is calculated by ’’Hard 
Thermal Loop” approximation as follows [42] : 


dE_ 

~dL 


8 a 2 s T 2 
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where rrijj = ^-(3 + n,)T 2 is the Debye mass. Finally the total energy loss due to both 
collision and radiation precesses is [43] : 
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where £ = 2 9 3t . Using equation ( 8 ) — (11), we will calculate the drag and diffusion 
coefficients for HQ energy loss during the interaction with QGP medium. In the next 
section we will find the solution of the Fokker-plank equation, eq( 6 ). 


III. The Iterative Laplace Transform Method 

In this section we want to introduce a new method for solving the FP equation, namely 
the ’’Iterative Laplace Transform Method” (ILTM) [44]. Many problems in physics and 
engineering can be successfully modelled by fractional differential equations (FDE). The 
FP equation indeed, is an important example of FDEs. Standard methods like finite- 
difference method are not applicable for such differential equations with non-zero values 
on the boundaries of initial conditions. Hence should we be looking for an effective way 
to solve FDEs. This new method was proposed by Daftardar-Gejji and Jafari to seek 
numerical solutions of nonlinear functional equations. This method is called the Iterative 
Laplace Transform Method, which is a combination of two powerful methods, namely, 
the Laplace transform method and the Iterative method. The method gives numerical 
solutions in the form of convergent series with easily computable components. The most 
outstanding feature of this method is that it provides an analytical solution using the 
initial conditions only, without any discretization or restrictive assumptions. 

In general, the Fokker-Planck equation can be written as [44] : 


D?f = 


D g A(p, t, f) + D 2 /B(p , t, f) 


f(p,t) 


( 12 ) 


Here Df (.), Dp (.), D 2 ^(.) are the Caputo fractional derivative with respect to t and p. 
In our calculations we use the following definitions: 

I. The Caputo fractional derivative of function (p, t ) is defined as 
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D?f(p, t) = 1 rir-^r(p, rj)drj (13) 

r (m - a) J 

for m — 1 < a < ?n,m € N. 

II. The Laplace transform of /(f) is defined as: 

roo 

F(s) = if [/(f)] = / e~ st f(t)dt . (14) 

Jo 

III. Laplace transform of D^f(p,t) is given by 

m— 1 

■^[ D tf{p,t)] = x a ^f[f(p,t)] - E f ( ' k \p,0)x a ~ 1 ~ k ; m — 1 < a < m , (15) 

k =0 

where f^ k \p, 0) is the &ith derivative of f(p,t ) at f = 0. 

To illustrate how the iterative Laplace transform method works, the general space-time 
fractional partial differential equation is considered 

D?f = s/{f,D$f,Dff,...) m- 1 < a < m 

n—l</3<n m,n£N , (16) 

where (/, Dpf , Dff ,...) is a linear or nonlinear operator of /, Dpf , D'f /,.... 

With the initial condition 

f (k \p,0) = h k (p)] k = 0, l,...,m - 1 , (17) 

and f{p,t) will be determined later. 

First of all we take the Laplace transform of both sides of (16) 

m— 1 

x a Sf[f{p, f)] - E x“- 1 - fc / (fe) (P, 0) = if W, £>£/, Df f ,...)] . (18) 

/c=0 

After simplification, we have 

m—1 

i?[/(p,t)] = E *" 1_fc / (fc) (p,0) +x-“if[^(/,^/, Df /,...)]. (19) 

fc=o 

By inverse Laplace transforms of both sides of Eq. (19) we have: 

m—1 

f(p, t ) = i^[E * _1_fc / (fc) (p, 0)] + if-V^a, £>£/, Df /,...)]] , (20) 

fc=0 

and 

#(/, D%f, Dff ,...) = if-Hx^if [^(/, D%f, Dff ,...)]] . (21) 

The method gives numerical solutions in the form of convergent series 


5 



( 22 ) 


f(P,t) = fn > 

n=0 


and 


Yfn,D 2 /Yf- 


n =0 


n =0 


= »(/0,O»/„,Df/„) 


+E 

j=o 


E 

j=i 



(23) 


By substituting (22) and (23) into (20) we have 


r m—1 


n =0 


Yfn= ■Z’- 1 0 ) 


L fc =0 

+@(f 0 ,DPf 0 ,D 2 /f 0 ...) 

OO r , j j j 

+E ®[Yfk, D iYfk> D ¥ E/fc.- 

j=i L Vfe=o fc=o fc=o 

Yfk,D^ P Yfk,D 2 /Yfk, 


lt =o 


fc =0 


fc =0 


Then we set 


(24) 


r m—1 


./o = ^ 


-1 


Y*- l ~ k f (k Xp, o) 

L fc =0 

fi=&(fo,Dlfo,DWf 0 ,...) 


fm+i = &\ Y^ d pY /*> D f E /*> - 

' fc=0 /c=0 fc=0 

-m—1 m—1 m—1 


E E D f E ^ 1 


fc =0 


k=0 


k =0 


( 25 ) 


Therefore the solution of (16)with initial condition (17)is given by: 
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f(p,t) = fo(p,t) + h(p,t) + ... + f m (p, t),m = 1,2,... (26) 

If we set a = j3 = 1 in Equation (26) we will reach the FPE. The drag and diffusion 
coefficients are explicitly timely independent parameters. However we recalculate these 
parameters before every time steps of our simulation. In this situation we have: 

f(p,t + At) = /o + /i + ... + f m (27) 


where 


fo 

fi 


f(p,t ) 

At^A(p)-^-(pf) + D(p) 


d_ 

dp 


f=fo 


fm+l — 


(At) 


m+1 


m + l 


A(p)-^(pf) + d (p) 


d_ 

dp 


, m > 1 . 


(28) 


f=fn 


We solved the FPE using both finite difference method and ILTM with different types of 
initial conditions and compared their results. For initial conditions with non vanishing 
values at the borderer, the ILTM successfully solves the FPE while results of finite dif¬ 
ference method diverges and therefore it failes to find the solution. It may be noted that 
momentum distribution functions of HQs, generally find their maximum at (p = 0). 

We have taken momentum distribution functions of charm and bottom quarks at yfs = 
200GeV in Au-Au and also P-P collisions from ’’The Durham HepData Project” database 
[45] which are simulated data. ILTM has been performed up to m = 2, which means we 
calculate fo, fo and fo in (128[) . 


IV. Time evolution of HQ distribution functions in QGP 

The dissipation coefficients of the thermal bath and initial momentum distribution of HQs, 
are basic inputs required to solve the FPE. Time evolution of dynamical parameters (like 
temperature, viscosity and so on) should be taken from suitable models. In our calculation 
the time dependence of temperature is considered as [46] 

T(t) = r 0 1/3 T 0 /t 1/3 , (29) 

where tq and To are the initial time and initial temperature that the background of the 
partonic system has attained in local equilibrium. We have taken tq = 0.33 fm/c and 
To = 0.373 GeV in our simulations [46]. The viscosity effects have not considered for 
simplification of our problem. Here we considered a fixed value for coupling constant as 
a s = 0.3. It is clear that the background is not stationary, it expands and cools down with 
time, so for more accurate results one can consider a strong coupling, a s , running with 
HQ mass and/or the medium temperature. For calculating the time evolution of the HQ 
momentum distribution by solving the FPE using ILTM, we need three inputs: I) initial 
distribution function, II) drag coefficient A(p) which can be calculated by Inserting (fill] 
in the Equation ([7]) and III) the diffusion coefficient as D = A(p)TM. Finally, we find the 
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momentum distribution of HQs at transition temperature T c = 175 MeV. It is clear that 
more advanced relations and more suitable assumptions can produce better results which 
can be found in further works. 

Figure 1 demonstrates the drag coefficient for b and c quarks due to collision and 
radiation separately. This figure clearly shows that the energy loss for c quark is greater 
than that for b quark. Figure 2 presents the HQ momentum distributions at initial time 



Figure 1: The drag coefficient due to collision and radiation as functions of momentum for Bottom and 
Charm quarks in QGP. 

(r = 0.33 fm/c) and at final time, when the QGP temperature becomes T = 0.175 GeV. 
A general comparison of momentum distributions before and after the evolution indicates 
that the probability of finding quarks with lower momentums after the evolution is greater 
than the probability of finding low momentum quarks at initial time. It may be noted 
that momentum distributions can not be compared directly, because they describe quark 
distributions at different values of total energy. This figure also shows that energy loss 
due to suppression of c quark is greater than that for b quark. 


V. Calculating the nuclear suppression factor Raa{Pt) 

One of the most important experimental observables to quantify the nuclear medium is the 
depletion of high px particles (D and B mesons or non-photonic single e _ spectra resulting 
from semi leptonic decays of hadrons containing charm and bottom quarks) produced 
in Nucleus-Nucleus collisions with respect to those produced in proton-proton collisions 
which is formulated in nuclear suppression factor Raa■ The nuclear suppression factor 
is calculated by solving relativistic hydrodynamic equations for the background medium 
in QGP phase, simultaneously with the FP equation for the HQs.The properties of the 
background bath affect the energy loss of HQs through the drag and diffusion coefficients. 
The energy loss of HQs is sensitive to initial conditions and the model used to describe 
the dynamics of the system. 
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Figure 2: Time evolution of b and c quark momentum distributions propagating in QGP medium. Initial 
values are: r 0 = 0.33^ , To = 0.373GeP and a s = 0.3. 


Here the nuclear modification factor ( Raa ) hi relativistic heavy ion collisions is calcu¬ 
lated for demonstrating the power of ILTM applied to solve the FP equation. 

We now proceed with the presentation of our numerical simulation with the formalism 
described and the inputs mentioned in the previous section. Then we present the compari¬ 
son between the results obtained with solving FPE by ILTM approach by the experimental 
data. There are several free parameters in the problem, which should be initiated or cho¬ 
sen to find an acceptable fitting the results on the reference data. Equilibration time after 
the collision, equilibration temperature, contribution of collisional and radiation dissipa¬ 
tions on drag and diffusion coefficients, effective light quark flavors iif , the charm and the 
bottom mass, hadron multiplicity at mid-rapidity ^ and final temperature are some of 
needed parameters. 

To make a realistic connection between simulation results and the experimental data, 
we have to implement the hadronization mechanisms at the quark-hadron transition tem¬ 
perature in order to find the non-photonic single electron spectra originating from the 
decays of heavy flavoured mesons (D and B). The solution of the FP equation (described 
in the previous section) have been convoluted with the fragmentation functions of the 
heavy quarks to obtain the Pt distribution of the heavy mesons. We have used the fol¬ 
lowing heavy quark fragmentation [561 [57] : 


f(z) « 


z\ z 



(30) 


We have assumed that the electron spectrum produced through semileptonic decay of D 
and B-mesons is proportional to the heavy meson spectra. Also we have not considered 
the expansion of the QGP bath. These assumptions create some expected deviation from 
the experimental data. 
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Similar simulations also have been performed to find the electron spectrum created in 
the P-P collisions. For this purpose HQ distribution functions of proton at y/s = 200 GeV 
have been taken as initial conditions to solve FPE again. The ratio of these two quantities 
is proportional to the nuclear suppression factor, Raa measured in experiments: 

djye Au—Au 

Raa(Pt) oc ^ p_p (31) 

dPTdr) 

A small value of Raa indicates a strong suppression and therefore, large energy loss of 
heavy quarks in the medium. It is clear that this ratio is equal to one in the absence of 
any medium. 

Figure 3 presents contribution of collision and radiation energy loss in suppression 
factor Raa separately. Displayed experimental data has been obtained from the PHENIX 
collaborations for Au + Au collisions at y/sjyN = 200 GeV. 



P T GeV/c 


Figure 3: Nuclear suppression factor Raa due to collisional and radiation energy loss separately as 
functions of Pt . Used parameters are: to = 0.33/m/e, To = 0.375GeV, a s = 0.3, T c = 0.175GeU. 

This figure shows that suppression at higher values of Pt can not be explained by 
collisional energy loss. In other words, contribution of radiation in energy loss at higher 
Pt is more important than collisional energy loss only. Figure 4 demonstrates our best 
simulation result by considering both collision and radiation energy loss in calculating the 
drag and diffusion coefficient to solve the FPE. The drag coefficient has been considered as 
Atotai = A coUison + 1.7 A radiation . This figure shows a very good agreement between results 
and experimental data at least for Pt < 5 GeV. It may be noted that the final temperature 
has been taken T c = 0.16 GeV and an offset 0.05 has been added too. As mentioned before 
we have not considered some details of the real situation in our simulations. 
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Figure 4: Nuclear suppression factor Raa as function of Pt ■ Energy loss due to collision and radiation 
have been considered in calculating the drag coefficient. 


IV. Conclusions and remarks 

The ’’Iterative Laplace Transform Method” has been introduced as an effective method 
to solve the Fokker-Planck equation for initial conditions with non-zero values at the 
boundaries. A numerical algorithm for solving FPE is presented in this paper. Our 
calculations show that this method is able to solve time evolution of distribution function 
of HQs successfully. To verify the ability of ILTM and demonstrating its application, we 
calculated the nuclear suppression factor Raa■ It is shown that there is a good agreement 
between simulation results and reported experimental data. 

Therefore many problems can be investigated using the ILTM. Effects of a running 
coupling constant, HQ mass, dead cone effect and LPM effects, QGP viscosity, different 
equations of state as models of QGP are subjects which can be studied in further works. 
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